#include "base.h"

int main(int argc,char **argv)
{


//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@///

	  //////////////////
	  // Input   ///////
	  //////////////////

	  int Nr=2;
	  char kappa_sea[1024];
	  int conf_ini, conf_end, jump,Nskip;
	  int skipList[50];
	  int Nconf, Nex;
	  int Nmass;
	  char dir_input[1024];
	  char dir_output[1024];
	  char P_to_use[1024];
	  int  set_p[4];
	  int iprop;

	  int twisted_basis;

	  int beta100;
          char gauge_fixing[1024];
          char gauge_action[1024];
          double Uplaq;
          int NOa2;
          double beta;

          char momentum_order[1024];
          char gamma[1024];

	  ////// Parsing command line

	 int iskip=0;
    	  for(int i=1; i<argc; i++) {
	    if(strncmp(argv[i],"-kappa_sea",10)==0)  sscanf(argv[i],"-kappa_sea=%s",kappa_sea);
	    if(strncmp(argv[i],"-conf_ini",9)==0)  sscanf(argv[i],"-conf_ini=%d",&conf_ini);
            if(strncmp(argv[i],"-jump",5)==0)   sscanf(argv[i],"-jump=%d",&jump);
	    if(strncmp(argv[i],"-Nskip",6)==0)  sscanf(argv[i],"-Nskip=%d",&Nskip);
	    if(strncmp(argv[i],"-skip",5)==0)  { sscanf(argv[i],"-skip=%d",&skipList[iskip]); iskip++; }
	    if(strncmp(argv[i],"-Nclust",7)==0)  sscanf(argv[i],"-Nclust=%d",&Nclust);
	    if(strncmp(argv[i],"-Nex",4)==0)  sscanf(argv[i],"-Nex=%d",&Nex);

	    if(strncmp(argv[i],"-Nmass",6)==0)  sscanf(argv[i],"-Nmass=%d",&Nmass);
            if(strncmp(argv[i],"-dir_input",10)==0)  sscanf(argv[i],"-dir_input=%s",dir_input);
            if(strncmp(argv[i],"-dir_output",11)==0)  sscanf(argv[i],"-dir_output=%s",dir_output);
	    if(strncmp(argv[i],"-T",2)==0)  sscanf(argv[i],"-T=%d",&T); 
            if(strncmp(argv[i],"-L",2)==0)  sscanf(argv[i],"-L=%d",&L);
            if(strncmp(argv[i],"-beta100",8)==0)  sscanf(argv[i],"-beta100=%d",&beta100);
	    if(strncmp(argv[i],"-P_to_use",9)==0)  sscanf(argv[i],"-P_to_use=%s",P_to_use); 
	    if(strncmp(argv[i],"-iprop",6)==0)  sscanf(argv[i],"-iprop=%d",&iprop);
	    if(strncmp(argv[i],"-set_p",6)==0)  sscanf(argv[i],"-set_p=%d,%d,%d,%d",&set_p[3], &set_p[2], &set_p[1], &set_p[0]);

            if(strncmp(argv[i],"-beta",5)==0)  sscanf(argv[i],"-beta=%lf",&beta);
            if(strncmp(argv[i],"-Uplaq",6)==0)  sscanf(argv[i],"-Uplaq=%lf",&Uplaq);
            if(strncmp(argv[i],"-NOa2",5)==0)  sscanf(argv[i],"-NOa2=%d",&NOa2);
            if(strncmp(argv[i],"-gauge_fixing",13)==0)  sscanf(argv[i],"-gauge_fixing=%s",gauge_fixing);
            if(strncmp(argv[i],"-gauge_action",13)==0)  sscanf(argv[i],"-gauge_action=%s",gauge_action);

             if(strncmp(argv[i],"-momentum_order",15)==0)  sscanf(argv[i],"-momentum_order=%s",momentum_order);
             if(strncmp(argv[i],"-gamma",6)==0)  sscanf(argv[i],"-gamma=%s",gamma);

             if(strncmp(argv[i],"-twisted_basis",13)==0)  sscanf(argv[i],"-twisted_basis=%d",&twisted_basis);

    	  }
	  Nconf=Nclust*Nex;
	  conf_end=conf_ini+jump*(Nconf-1)+Nskip*jump;
	  int Nmass_comb=Nmass;
       	  for (int i=1; i<Nmass; i++)  Nmass_comb=Nmass_comb+(Nmass-i);
          init_base_gamma(gamma);

          printf("==============================\n");
          printf("conf_ini=%d, conf_end=%d jump=%d\nNconf=%d, Nclust=%d, Nex=%d\n",conf_ini, conf_end,jump,Nconf,Nclust,Nex);
          printf("Nmass=%d\n",Nmass);
	  printf("Nmass_comb=%d\n",Nmass_comb);
          printf("L=%d T=%d\n",L,T);
          printf("Nr=%d \n",Nr);
	  printf("set=(%d,%d,%d,%d)\n",set_p[0], set_p[1], set_p[2], set_p[3]);
	  if (twisted_basis) printf(" Propagators in base twisted\n");
	  else  printf(" Propagators in physical base \n");
          printf("dir_input=%s \n",dir_input);
          printf("dir_output=%s \n",dir_output);
	  if (NOa2>1) {
                        printf("O(a2g2) Corrections beta=%0.3f -Uplaq=%0.5f %s %s\n", beta, Uplaq, gauge_fixing, gauge_action);
         }

          printf("==============================\n");
	
   /////////////////////////
  //// Momenta set
  //////////////////////////////
  char momentafile[1024]; sprintf(momentafile,"../INPUT/mom_%d.dat",beta100);
  read_momenta_data(momentafile);
  double *nn=(double*)malloc(sizeof(double)*(Nmom));
  compute_nn(momentum_order,nn);
  Momentum_Magro(momentum_order,set_p); 
  // equivalent momenta
  Nmom_eq=Number_of_Equivalent_Momenta(SinP2);
  printf("Nmom_eq %d \n",Nmom_eq);

  Average_over_Equivalent_Momenta(SinP2,SinP2,SinP2_eq);


  double *P2tilde=(double*)malloc(sizeof(double)*(Nmom));
  double **Ptilde=(double**)malloc(sizeof(double*)*(Nmom));
  for (int imom=0 ; imom<Nmom; imom++) Ptilde[imom]=(double*)malloc(sizeof(double)*(4));

  for (int imom=0 ; imom<Nmom; imom++) P2tilde[imom]=SinP2[imom];
  for (int imom=0 ; imom<Nmom; imom++)  for (int idir=0; idir<4; idir++) Ptilde[imom][idir]=SinP[imom][idir];
  if(strcmp(P_to_use,"P")==0)
  {
                 printf("Using P insteadd SinP\n");
                  for (int imom=0 ; imom<Nmom; imom++) P2tilde[imom]=P2[imom];
                  for (int imom=0 ; imom<Nmom; imom++)  for (int idir=0; idir<4; idir++) Ptilde[imom][idir]=P[imom][idir];

  }
  

	  //////////////////////////////////
          // IO files
          ///////////////////////////////

          char s0ff_base_input[1024];
          char s0ft_input[1024];
          char z_base_output[1024];

	  char file[1024];

	  char gammaprint[8];
	  char Oa2print[8];
          char dir_qav[1024];
          sprintf(dir_qav,"%s/QAV",dir_output);
 
          sprintf(s0ff_base_input,"%s/s%dft",dir_input,iprop);

          sprintf(z_base_output,"%s/Z_ALLP/Z",dir_output);


//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@//

           /////////////////////////////
	   // Read Fourier transform
	   //////////////////////////////
	
          //def
	  colorcolorspinspin ****s0FT_r=malloc(sizeof(colorcolorspinspin***)*(Nconf));
          for (int iconf=0; iconf<Nconf; iconf++) s0FT_r[iconf]=malloc(sizeof(colorcolorspinspin**)*(Nmass));
          for (int iconf=0; iconf<Nconf; iconf++)  for(int imass=0; imass<Nmass; imass++) s0FT_r[iconf][imass]=malloc(sizeof(colorcolorspinspin*)*(Nmom)); 
	  for (int iconf=0; iconf<Nconf; iconf++)  for(int imass=0; imass<Nmass; imass++) for (int imom=0; imom<Nmom; imom++) s0FT_r[iconf][imass][imom]=malloc(sizeof(colorcolorspinspin)*(Nr));

	  int conf_to_read;
	  FILE *fin=NULL;


           printf("====================\n");
           printf("READING...\n");

           for (int r=0; r<Nr; r++)
	   {
	   int conf_skiped_up_to_now=0;

           for (int iconf=0; iconf<Nconf+Nskip; iconf++)
 	   {
                   conf_to_read=conf_ini+iconf*jump;
		   int skip_this_conf=0;
	         for (int iskip=0; iskip<Nskip; iskip++)
		              if (conf_to_read==skipList[iskip]){
	                      skip_this_conf=1;
                 	      conf_skiped_up_to_now++;
                               }
                 if( skip_this_conf==1 )  printf("Skiping conf %d\n",conf_to_read); 
		 else{
		   int iconf_real=iconf-conf_skiped_up_to_now;
                   //printf("  ... conf %d\n",conf_to_read);
                   sprintf(s0ft_input,"%s%d_%s_%04d.out",s0ff_base_input,r,kappa_sea,conf_to_read);
                   fin=fopen(s0ft_input,"rb+");
                   if(fin == (FILE *)NULL) {
                        //fprintf(stderr,"Unable to open file %s for reading\n",s0ft_input);
                        sprintf(s0ft_input,"%s%d_%04d.out",s0ff_base_input,r,conf_to_read);
                        fin=fopen(s0ft_input,"rb+");
                        if(fin == (FILE *)NULL)
                        {
                                //fprintf(stderr,"Unable to open file %s for reading\n",s0ft_input);
                                sprintf(s0ft_input,"%s%d_%s_%03d.out",s0ff_base_input,r,kappa_sea,conf_to_read);
                                fin=fopen(s0ft_input,"rb+");
                                if(fin == (FILE *)NULL)
                                {
                                        fprintf(stderr, "Unable to open file %s%d_%s_%04d.out  or %s%d_%04d.out or %s%d_%s_%03d.out",s0ff_base_input,r,kappa_sea,conf_to_read,s0ff_base_input,r,conf_to_read,s0ff_base_input,r,kappa_sea,conf_to_read );
                                        return EXIT_FAILURE;

                                }
                        }
                   }
                        printf("Reading %s\n",s0ft_input);

	          for(int imass=0; imass<Nmass; imass++)
                   {
                  // printf("im=%d\n",imass);
		   for (int imom=0; imom<Nmom; imom++)
			for (int ic_ext=0; ic_ext<3; ic_ext++) for (int is_ext=0; is_ext<4; is_ext++) for (int ic_int=0; ic_int<3; ic_int++) for (int is_int=0; is_int<4; is_int++){
			fread(&s0FT_r[iconf_real][imass][imom][r][ic_int][ic_ext][is_int][is_ext],sizeof(complex),1,fin);
			//printf("%8.10f +i %8.10f\n", creal(s0FT_r[iconf][imass][imom][r][ic_int][ic_ext][is_int][is_ext]),cimag(s0FT_r[iconf][imass][imom][r][ic_int][ic_ext][is_int][is_ext]));
			}
		   }
		   fclose(fin);
		}
	       }
	   }
	   printf("====================\n\n\n");


	    ////////////////////////////////////
            /////         CLUSTERIZE      //////
	    ////////////////////////////////////
	    printf("====================\n");
            printf("CLUSTERIZE\n");

	   complex   aux_conf[Nconf], aux_clust[Nclust+1];

          //def
          colorcolorspinspin ****s0FT=malloc(sizeof(colorcolorspinspin***)*(Nclust+1));
          for (int iclust=0; iclust<Nclust+1; iclust++) s0FT[iclust]=malloc(sizeof(colorcolorspinspin**)*(Nmass));
          for (int iclust=0; iclust<Nclust+1; iclust++)  for(int imass=0; imass<Nmass; imass++) s0FT[iclust][imass]=malloc(sizeof(colorcolorspinspin*)*(Nmom));
	  for (int iclust=0; iclust<Nclust+1; iclust++)  for(int imass=0; imass<Nmass; imass++) for (int imom=0; imom<Nmom; imom++) s0FT[iclust][imass][imom]=malloc(sizeof(colorcolorspinspin)*(Nr));

	  for (int r=0; r<Nr; r++)
	    for(int imass=0; imass<Nmass; imass++)  for (int imom=0; imom<Nmom; imom++)
            {
		for (int ic_ext=0; ic_ext<3; ic_ext++) for (int is_ext=0; is_ext<4; is_ext++) for (int ic_int=0; ic_int<3; ic_int++) for (int is_int=0; is_int<4; is_int++)
		{	
	               for(int iconf=0;iconf<Nconf;iconf++) aux_conf[iconf]=s0FT_r[iconf][imass][imom][r][ic_int][ic_ext][is_int][is_ext];
                       conf2clust(aux_conf, aux_clust,Nconf,Nclust,Nex);
                       for(int iclust=0;iclust<Nclust+1;iclust++) s0FT[iclust][imass][imom][r][ic_int][ic_ext][is_int][is_ext]=aux_clust[iclust];
		}
            }
	   
            printf("====================\n\n\n");



//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@//

	/////////////////////////////////////////////////////
        // FORM FACTOR OF THE FOURIER-TRANSFORMED PROPAGATOR
	//////////////////////////////////////////////////////

        printf("====================\n");
        printf("FORM FACTOR\n");

	//def
	complex *****s0_ff=malloc(sizeof(complex****)*(Nclust+1));
        for (int iclust=0; iclust<Nclust+1; iclust++) s0_ff[iclust]=malloc(sizeof(complex***)*(Nmass));
        for (int iclust=0; iclust<Nclust+1; iclust++) for(int imass=0; imass<Nmass; imass++) s0_ff[iclust][imass]=malloc(sizeof(complex**)*(Nmom));
	 for (int iclust=0; iclust<Nclust+1; iclust++) for(int imass=0; imass<Nmass; imass++) for (int imom=0; imom<Nmom; imom++) s0_ff[iclust][imass][imom]=malloc(sizeof(complex*)*(Nr));
	for (int iclust=0; iclust<Nclust+1; iclust++) for(int imass=0; imass<Nmass; imass++) for (int imom=0; imom<Nmom; imom++) for (int r=0; r<Nr; r++) s0_ff[iclust][imass][imom][r]=malloc(sizeof(complex)*(3));
	

	 for (int imom=0; imom<Nmom; imom++)
	 {
		        for (int iclust=0; iclust<Nclust+1; iclust++) for (int imass=0; imass<Nmass; imass++) for (int r=0; r<Nr; r++){
			compute_s0tw_ff(s0FT[iclust][imass][imom][r],s0_ff[iclust][imass][imom][r],Ptilde[imom],P2tilde[imom],nn[imom],r);
			}
	  }


        printf("====================\n\n\n");
        

	////////////////////////////////
	// INVERTING		////////
	////////////////////////////////

        printf("====================\n");
	printf("INVERTING PROPAGATOR\n");

        colorcolorspinspin ****s0FT_inv=malloc(sizeof(colorcolorspinspin***)*(Nclust+1));
        for (int iclust=0; iclust<Nclust+1; iclust++) s0FT_inv[iclust]=malloc(sizeof(colorcolorspinspin**)*(Nmass));
        for (int iclust=0; iclust<Nclust+1; iclust++)  for(int imass=0; imass<Nmass; imass++) s0FT_inv[iclust][imass]=malloc(sizeof(colorcolorspinspin*)*(Nmom));
	for (int iclust=0; iclust<Nclust+1; iclust++)  for(int imass=0; imass<Nmass; imass++) for (int imom=0; imom<Nmom; imom++) s0FT_inv[iclust][imass][imom]=malloc(sizeof(colorcolorspinspin)*(Nr));



	complex **mdum=(complex**)malloc(sizeof(complex)*(12));
	for(int i=0; i<12; i++) mdum[i]=(complex*)malloc(sizeof(complex)*(12));

        complex **idum=(complex**)malloc(sizeof(complex)*(12));
        for(int i=0; i<12; i++) idum[i]=(complex*)malloc(sizeof(complex)*(12));

	colorcolorspinspin s0_dump;

	int ind_ext, ind_int;
	int iflag=0; 

	for (int r=0; r<Nr; r++)
        for (int imom=0; imom<Nmom; imom++)
        {
	    for (int iclust=0; iclust<Nclust+1; iclust++) for (int imass=0; imass<Nmass; imass++)
            {
                        for (int ic_ext=0; ic_ext<3; ic_ext++) for (int is_ext=0; is_ext<4; is_ext++) for (int ic_int=0; ic_int<3; ic_int++) for (int is_int=0; is_int<4; is_int++)
                         {
				ind_ext=is_ext+4*(ic_ext);
				ind_int=is_int+4*(ic_int);
				mdum[ind_int][ind_ext]=s0FT[iclust][imass][imom][r][ic_int][ic_ext][is_int][is_ext];			

			}
             		inverse(mdum,12,idum,iflag);
			if (iflag==1) {break; printf("Light quark propagator is singular\n");}
			for (int ic_ext=0; ic_ext<3; ic_ext++) for (int is_ext=0; is_ext<4; is_ext++) for (int ic_int=0; ic_int<3; ic_int++) for (int is_int=0; is_int<4; is_int++)
                         {
                                ind_ext=is_ext+4*(ic_ext);
                                ind_int=is_int+4*(ic_int);
				s0FT_inv[iclust][imass][imom][r][ic_int][ic_ext][is_int][is_ext]=idum[ind_int][ind_ext];
			}
			/////////////////////////
			// Inversion Test
			///////////////////////
			for (int ic_ext=0; ic_ext<3; ic_ext++) for (int is_ext=0; is_ext<4; is_ext++) for (int ic_int=0; ic_int<3; ic_int++) for (int is_int=0; is_int<4; is_int++)
			{
				s0_dump[ic_int][ic_ext][is_int][is_ext]=0.;
				for (int ic_dump=0; ic_dump<3; ic_dump++) for (int is_dump=0; is_dump<4; is_dump++)
						s0_dump[ic_int][ic_ext][is_int][is_ext]=s0_dump[ic_int][ic_ext][is_int][is_ext]+s0FT[iclust][imass][imom][r][ic_int][ic_dump][is_int][is_dump]*s0FT_inv[iclust][imass][imom][r][ic_dump][ic_ext][is_dump][is_ext];
			}
			for (int ic_ext=0; ic_ext<3; ic_ext++) for (int is_ext=0; is_ext<4; is_ext++) for (int ic_int=0; ic_int<3; ic_int++) for (int is_int=0; is_int<4; is_int++)
			{
				if (ic_ext != ic_int || is_ext != is_int) 
				{
					if (cabs(s0_dump[ic_int][ic_ext][is_int][is_ext])> 0.1) 
						printf("Inversion failed for imom=%d (s0*s0inv)_%d%d%d%d=%8.5f+I%8.5f \n",imom,ic_int,ic_ext,is_int,is_ext,creal(s0_dump[ic_int][ic_ext][is_int][is_ext]),cimag(s0_dump[ic_int][ic_ext][is_int][is_ext]));
				}
				else
                                 {
				       if (cabs(s0_dump[ic_int][ic_ext][is_int][is_ext])-1> 0.1) 
						printf("Inversion failed  for imom=%d (s0*s0inv)_%d%d%d%d=%8.5f+I%8.5f \n",imom,ic_int,ic_ext,is_int,is_ext,creal(s0_dump[ic_int][ic_ext][is_int][is_ext]),cimag(s0_dump[ic_int][ic_ext][is_int][is_ext]));
				}
				
			}
                        ////////////////////////////////////
                        //  end inversion test
                        ////////////////////////////////////

	     }
        }
	printf("====================\n\n\n");

       /////////////////////////////////////////////////////
        // FORM FACTOR OF THE INVERTED FOURIER-TRANSFORMED PROPAGATOR
        //////////////////////////////////////////////////////

	printf("====================\n");
        printf("FORM FACTOR OF THE INVERTED PROPAGATOR\n");
        complex *****s0_inv_ff=malloc(sizeof(complex****)*(Nclust+1));
        for (int iclust=0; iclust<Nclust+1; iclust++) s0_inv_ff[iclust]=malloc(sizeof(complex***)*(Nmass));
        for (int iclust=0; iclust<Nclust+1; iclust++) for(int imass=0; imass<Nmass; imass++) s0_inv_ff[iclust][imass]=malloc(sizeof(complex**)*(Nmom));
	for (int iclust=0; iclust<Nclust+1; iclust++) for(int imass=0; imass<Nmass; imass++) for (int imom=0; imom<Nmom; imom++) s0_inv_ff[iclust][imass][imom]=malloc(sizeof(complex*)*(Nr));
        for (int iclust=0; iclust<Nclust+1; iclust++) for(int imass=0; imass<Nmass; imass++) for (int imom=0; imom<Nmom; imom++) for (int r=0; r<Nr; r++) s0_inv_ff[iclust][imass][imom][r]=malloc(sizeof(complex)*(3));


       for (int r=0; r<Nr; r++)
       for (int imom=0; imom<Nmom; imom++)
       {
		for (int iclust=0; iclust<Nclust+1; iclust++) for (int imass=0; imass<Nmass; imass++)
                compute_s0tw_inv_ff(s0FT_inv[iclust][imass][imom][r],s0_inv_ff[iclust][imass][imom][r],Ptilde[imom],P2tilde[imom],nn[imom],r);
       }





	printf("====================\n\n\n");
  
        /////////////////////////////////////////////////////
        // Z_q=-i/12 * Tr(sum_mu gamma_mu*SinP(mu)*S^-1(p))/P2tilde

        //////////////////////////////////////////////////////

        printf("====================\n");
        printf("Z_q\n");

        complex ****Zq=malloc(sizeof(complex***)*(Nclust+1));
        for (int iclust=0; iclust<Nclust+1; iclust++) Zq[iclust]=malloc(sizeof(complex**)*(Nmass));
        for (int iclust=0; iclust<Nclust+1; iclust++) for(int imass=0; imass<Nmass; imass++) Zq[iclust][imass]=malloc(sizeof(complex*)*(Nmom));
	for (int iclust=0; iclust<Nclust+1; iclust++) for(int imass=0; imass<Nmass; imass++) for (int imom=0; imom<Nmom; imom++) Zq[iclust][imass][imom]=malloc(sizeof(complex)*(Nr));


	for (int r=0; r<Nr; r++)
        for (int iclust=0; iclust<Nclust+1; iclust++) for (int imass=0; imass<Nmass; imass++) for (int imom=0; imom<Nmom; imom++)
	{
			Zq[iclust][imass][imom][r]=s0_inv_ff[iclust][imass][imom][r][0];
	}

       printf("====================\n\n\n");

       /////////////////////////////////////////////////////
       // Lambda_O'= Tr [ G_amp_O' * P_O'] ; O' = e^(-ig5 pi/4) O e(ig5 pi/4) 
       // Zq^-1 Z_O * Lambda_O' = 1 
       // The appeareance of O' is a trivial consequence of the fact that RCs,  are named after the (here unphysical) twisted quark basis while the operator OΓ = uΓd is expressed in the physical quark basis.
       // 
       // Lambda= Tr [ G_amp_O' * P_O' ] ; 
       // Twisted Mass combinations : r1=-r2 --> u,d
       // G_amp_O' = S_u^-1 G_O' S_d^-1 = S_u^-1 G_O' g5 S_u^\dagger g5
       // G_O'=<u(x) O' \bar d(x) > =S_u O' S_d^back=S_u O' g5 S_u^\dagger g5
       // Using TM fermions ZA <-> ZV
       // BUT if instead O' I use O both effects compensate: A'=-iV; V'=iA; P'=P; S'=S; T'=T
	// and i factors compensate with P_O=conj(O)

	//////////////////////////////////////////////////
	// So I use Lambda_O and I then I can read directly ZA and ZV
	/////////////////////////////////////////////////

       // P_O projector onto tree-level
       // if O=I,g5,g_mu,g5*g_mu... P_O=conj(O)
       //////////////////////////////////////////////////////

        printf("====================\n");
        printf("LAMBDA AND Z_O\n");

      complex *****Lambda=malloc(sizeof(complex****)*(5));
      for (int igamma=0; igamma<5 ; igamma++) Lambda[igamma]=malloc(sizeof(complex***)*(Nclust+1));
      for (int igamma=0; igamma<5 ; igamma++) for (int iclust=0; iclust<Nclust+1; iclust++) Lambda[igamma][iclust]=malloc(sizeof(complex**)*(Nmass_comb));
      for (int igamma=0; igamma<5 ; igamma++) for (int iclust=0; iclust<Nclust+1; iclust++) for(int imass=0; imass<Nmass_comb; imass++) Lambda[igamma][iclust][imass]=malloc(sizeof(complex*)*(Nmom));
      for (int igamma=0; igamma<5 ; igamma++) for (int iclust=0; iclust<Nclust+1; iclust++) for(int imass=0; imass<Nmass_comb; imass++) for (int imom=0; imom<Nmom; imom++)  Lambda[igamma][iclust][imass][imom]=malloc(sizeof(complex)*(Nr));

	complex Lambda_tmp[16];

      // auxiliares (solo por claridad de las formulas)

      colorcolorspinspin **G_r=malloc(sizeof(colorcolorspinspin*)*(16));
      for (int igamma=0; igamma<16 ; igamma++) G_r[igamma]=malloc(sizeof(colorcolorspinspin)*(Nconf));
      colorcolorspinspin **G=malloc(sizeof(colorcolorspinspin*)*(16));
      for (int igamma=0; igamma<16 ; igamma++) G[igamma]=malloc(sizeof(colorcolorspinspin)*(Nclust+1));
      colorcolorspinspin S_xy, S_yx, Sinv_xy, Sinv_yx;
      for (int r=0; r<Nr; r++)
      {
       printf("r=%d\n",r);
      for (int imom=0; imom<Nmom; imom++)
      {
	   printf("imom=%d\n",imom);
	   int imass=0;
	   for (int im2=0; im2<Nmass; im2++) for (int im1=im2; im1<Nmass; im1++)
	   {
		printf("m1=%d m2 =%d imass=%d\n",im1,im2, imass);
		//////////////////////////////////////////////
                // Compute NON amputed Green Function
		/////////////////////////////////////////////

		// On one single configuration the quark propagator is not transationally invariant. It is only recovered after averaging over the configurations
		// Then we compute G over the initial configurations 
		// and then we average
                // G(p)=\sum_n G_n(p)/N
//		printf("- Compute Non Amputed GF\n");
		for (int iconf=0; iconf<Nconf; iconf++) 
	        {
			// 1. Initialize
			for (int ic1=0; ic1<3 ; ic1++) for (int ic2=0; ic2<3 ; ic2++) 
                        for (int is1=0; is1<4 ; is1++) for (int is2=0; is2<4 ; is2++) 
			{
				S_xy[ic1][ic2][is1][is2]=s0FT_r[iconf][im1][imom][r][ic1][ic2][is1][is2];
				S_yx[ic1][ic2][is1][is2]=((base_gamma[5]).entr[is1])*conj(s0FT_r[iconf][im2][imom][r][ic2][ic1][(base_gamma[5]).pos[is2]][(base_gamma[5]).pos[is1]])*(base_gamma[5]).entr[is2];
			}
			// 2. Compute
          	      for (int igamma=0; igamma<16; igamma++) {
                       for (int ic1=0; ic1<3 ; ic1++) for (int ic2=0; ic2<3 ; ic2++)
                       for (int is1=0; is1<4 ; is1++) for (int is2=0; is2<4 ; is2++)
			{
				 G_r[igamma][iconf][ic1][ic2][is1][is2]=0.;
				 for (int ic3=0; ic3<3 ; ic3++) for (int is3=0; is3<4 ; is3++){
				 G_r[igamma][iconf][ic1][ic2][is1][is2]=G_r[igamma][iconf][ic1][ic2][is1][is2]+S_xy[ic1][ic3][is1][is3]*((base_gamma[igamma]).entr[is3])*S_yx[ic3][ic2][(base_gamma[igamma]).pos[is3]][is2];
			 	  }
			}
		    }
			
		}
		
		//3. Clusterize
//		printf("- Clusterize \n");
		for (int igamma=0; igamma<16; igamma++)
                for (int ic1=0; ic1<3 ; ic1++) for (int ic2=0; ic2<3 ; ic2++)
                for (int is1=0; is1<4 ; is1++) for (int is2=0; is2<4 ; is2++)
                {
                       for(int iconf=0;iconf<Nconf;iconf++) aux_conf[iconf]=G_r[igamma][iconf][ic1][ic2][is1][is2];
                       conf2clust(aux_conf, aux_clust,Nconf,Nclust,Nex);
                       for(int iclust=0;iclust<Nclust+1;iclust++) G[igamma][iclust][ic1][ic2][is1][is2]=aux_clust[iclust];
                }
		/////////////////////////////
                // Compute Lambda 
		///////////////////////////
//		printf( "- Compute Lambda \n");
		 for (int iclust=0; iclust<Nclust+1; iclust++)
                {

			// 1. Initialize
			for (int ic1=0; ic1<3 ; ic1++) for (int ic2=0; ic2<3 ; ic2++)
                        for (int is1=0; is1<4 ; is1++) for (int is2=0; is2<4 ; is2++)
                        {
				Sinv_xy[ic1][ic2][is1][is2]=s0FT_inv[iclust][im1][imom][r][ic1][ic2][is1][is2];
				Sinv_yx[ic1][ic2][is1][is2]=((base_gamma[5]).entr[is1])*conj(s0FT_inv[iclust][im2][imom][r][ic2][ic1][(base_gamma[5]).pos[is2]][(base_gamma[5]).pos[is1]])*(base_gamma[5]).entr[is2];
			}
			// 2. Compute
			for (int igamma=0; igamma<16; igamma++)
			{
			 Lambda_tmp[igamma]=0.;
			 for (int ic1=0; ic1<3 ; ic1++) for (int ic2=0; ic2<3 ; ic2++) for (int ic3=0; ic3<3 ; ic3++)
			 for (int is1=0; is1<4 ; is1++) for (int is2=0; is2<4 ; is2++) for (int is3=0; is3<4 ; is3++)
			 {
			  Lambda_tmp[igamma]=Lambda_tmp[igamma]+Sinv_xy[ic1][ic2][is1][is2]*G[igamma][iclust][ic2][ic3][is2][is3]*Sinv_yx[ic3][ic1][is3][(base_gamma[igamma]).pos[is1]]*conj((base_gamma[igamma]).entr[is1])/12.;
			 }
			}
			Lambda[0][iclust][imass][imom][r]=Lambda_tmp[0];
			Lambda[1][iclust][imass][imom][r]=(Lambda_tmp[1]+Lambda_tmp[2]+Lambda_tmp[3]+Lambda_tmp[4])/4.;
			Lambda[2][iclust][imass][imom][r]=Lambda_tmp[5];
			Lambda[3][iclust][imass][imom][r]=(Lambda_tmp[6]+Lambda_tmp[7]+Lambda_tmp[8]+Lambda_tmp[9]+Lambda_tmp[10]+Lambda_tmp[11])/6.;
			Lambda[4][iclust][imass][imom][r]=(Lambda_tmp[12]+Lambda_tmp[13]+Lambda_tmp[14]+Lambda_tmp[15])/4.;
		}
	    imass++;	
	    }	

      }
     }

     printf("====================\n");
     printf(" PROMEDIUM r0 r1   \n");

     double ****Zqr=malloc(sizeof(double***)*(Nclust+1));
     for (int iclust=0; iclust<Nclust+1; iclust++) Zqr[iclust]=malloc(sizeof(double**)*(Nmass));
     for (int iclust=0; iclust<Nclust+1; iclust++) for(int imass=0; imass<Nmass; imass++) Zqr[iclust][imass]=malloc(sizeof(double*)*(Nmom));
      for (int iclust=0; iclust<Nclust+1; iclust++) for(int imass=0; imass<Nmass; imass++) for (int imom=0; imom<Nmom; imom++) Zqr[iclust][imass][imom]=malloc(sizeof(double)*(NOa2));
     
	 double *****Lambdar=malloc(sizeof(double****)*(5));
      for (int igamma=0; igamma<5 ; igamma++) Lambdar[igamma]=malloc(sizeof(double***)*(Nclust+1));
      for (int igamma=0; igamma<5 ; igamma++) for (int iclust=0; iclust<Nclust+1; iclust++) Lambdar[igamma][iclust]=malloc(sizeof(double**)*(Nmass_comb));
      for (int igamma=0; igamma<5 ; igamma++) for (int iclust=0; iclust<Nclust+1; iclust++) for(int imass=0; imass<Nmass_comb; imass++) Lambdar[igamma][iclust][imass]=malloc(sizeof(double*)*(Nmom));
      for (int igamma=0; igamma<5 ; igamma++) for (int iclust=0; iclust<Nclust+1; iclust++) for(int imass=0; imass<Nmass_comb; imass++)  for (int imom=0; imom<Nmom; imom++) Lambdar[igamma][iclust][imass][imom]=malloc(sizeof(double)*(NOa2));

      double *****ZOr=malloc(sizeof(double****)*(5));
      for (int igamma=0; igamma<5 ; igamma++) ZOr[igamma]=malloc(sizeof(double***)*(Nclust+1));
      for (int igamma=0; igamma<5 ; igamma++) for (int iclust=0; iclust<Nclust+1; iclust++) ZOr[igamma][iclust]=malloc(sizeof(double**)*(Nmass_comb));
      for (int igamma=0; igamma<5 ; igamma++) for (int iclust=0; iclust<Nclust+1; iclust++) for(int imass=0; imass<Nmass_comb; imass++) ZOr[igamma][iclust][imass]=malloc(sizeof(double*)*(Nmom));
      for (int igamma=0; igamma<5 ; igamma++) for (int iclust=0; iclust<Nclust+1; iclust++) for(int imass=0; imass<Nmass_comb; imass++) for (int imom=0; imom<Nmom; imom++) ZOr[igamma][iclust][imass][imom]=malloc(sizeof(double)*(NOa2));


   for (int iclust=0; iclust<Nclust+1; iclust++) for (int imass=0; imass<Nmass; imass++) for (int imom=0; imom<Nmom; imom++)
        {
                Zqr[iclust][imass][imom][0]=(Zq[iclust][imass][imom][0]+Zq[iclust][imass][imom][1])/2.;
        }

     for (int igamma=0; igamma<5 ; igamma++)for (int iclust=0; iclust<Nclust+1; iclust++)  for (int imom=0; imom<Nmom; imom++)
     {
        int imass=0;
         for (int im2=0; im2<Nmass; im2++) for (int im1=im2; im1<Nmass; im1++)
           {
                Lambdar[igamma][iclust][imass][imom][0]=(Lambda[igamma][iclust][imass][imom][0]+Lambda[igamma][iclust][imass][imom][1])/2.;
                imass++;
           }
     }


     printf("====================\n\n\n");

     printf("====================================================\n");
     printf("O(a2g2) corrections \n ");


     //////////////// substracting terms     
     double ***CorOa2_Zq=malloc(sizeof(double**)*(NOa2));
     for (int iOa2=0; iOa2<NOa2; iOa2++) CorOa2_Zq[iOa2]=malloc(sizeof(double*)*(Nmom));
     for (int iOa2=0; iOa2<NOa2; iOa2++)  for (int imom=0; imom<Nmom; imom++) CorOa2_Zq[iOa2][imom]=malloc(sizeof(double)*(1));
     double ***CorOa2_Lambda=malloc(sizeof(double**)*(NOa2));
     for (int iOa2=0; iOa2<NOa2; iOa2++) CorOa2_Lambda[iOa2]=malloc(sizeof(double*)*(Nmom));
     for (int iOa2=0; iOa2<NOa2; iOa2++)  for (int imom=0; imom<Nmom; imom++) CorOa2_Lambda[iOa2][imom]=malloc(sizeof(double)*(5));

     //iOa2=0 -> uncorrected
     for (int imom=0; imom<Nmom; imom++){
     CorOa2_Zq[0][imom][0]=0.;
     for (int igamma=0; igamma<5 ; igamma++) CorOa2_Lambda[0][imom][igamma]=0.;
     }
     //iOa2=1 corrected
     //iOa2=2 corrected boosted
     double g2;
     for (int iOa2=1; iOa2<NOa2; iOa2++)
     {
        for (int imom=0; imom<Nmom; imom++)
        {
          if(iOa2==1) g2=6.0/beta;
          if(iOa2==2) g2=(6.0/beta)/Uplaq;
          correction_Oa2g2(CorOa2_Zq[iOa2][imom],CorOa2_Lambda[iOa2][imom],g2,SinP2[imom],SinP4[imom],nn[imom],gauge_action,gauge_fixing,twisted_basis);
        }

     }
      ////////////////// Zq corrected
      for (int iclust=0; iclust<Nclust+1; iclust++) for (int imom=0; imom<Nmom; imom++)
      for (int imass=0; imass<Nmass; imass++)
      {
         for (int iOa2=1; iOa2<NOa2; iOa2++) Zqr[iclust][imass][imom][iOa2]=Zqr[iclust][imass][imom][0]-CorOa2_Zq[iOa2][imom][0];
      }

      ///////////////// Lambda corrected
      for (int igamma=0; igamma<5 ; igamma++)
      for (int iOa2=1; iOa2<NOa2; iOa2++)
      for (int imass=0; imass<Nmass_comb; imass++)
        for (int iclust=0; iclust<Nclust+1; iclust++) for (int imom=0; imom<Nmom; imom++){
          Lambdar[igamma][iclust][imass][imom][iOa2]=Lambdar[igamma][iclust][imass][imom][0]-CorOa2_Lambda[iOa2][imom][igamma];
        }

    ///////////////// Z corrected

     for (int iOa2=0; iOa2<NOa2; iOa2++)
     for (int igamma=0; igamma<5 ; igamma++)for (int iclust=0; iclust<Nclust+1; iclust++)  for (int imom=0; imom<Nmom; imom++)
     {
        int imass=0;
         for (int im2=0; im2<Nmass; im2++) for (int im1=im2; im1<Nmass; im1++)
           {
                ZOr[igamma][iclust][imass][imom][iOa2]=sqrt(Zqr[iclust][im1][imom][iOa2]*Zqr[iclust][im2][imom][iOa2])/Lambdar[igamma][iclust][imass][imom][iOa2];
                imass++;
           }
     }





     printf("====================\n\n\n");

     //////////////////////////////////////////////
     //     EQUIVALENT MOMENTA
     ///////////////////////////////////////////////
     printf("==================================\n");
     printf(" PROMEDIUM OVER EQUIVALENT SINP2 \n");


        complex *****s0_ff_eq=malloc(sizeof(complex****)*(Nclust+1));
        for (int iclust=0; iclust<Nclust+1; iclust++) s0_ff_eq[iclust]=malloc(sizeof(complex***)*(Nmass));
        for (int iclust=0; iclust<Nclust+1; iclust++) for(int imass=0; imass<Nmass; imass++) s0_ff_eq[iclust][imass]=malloc(sizeof(complex**)*(Nmom_eq));
         for (int iclust=0; iclust<Nclust+1; iclust++) for(int imass=0; imass<Nmass; imass++) for (int imom=0; imom<Nmom_eq; imom++) s0_ff_eq[iclust][imass][imom]=malloc(sizeof(complex*)*(Nr));
        for (int iclust=0; iclust<Nclust+1; iclust++) for(int imass=0; imass<Nmass; imass++) for (int imom=0; imom<Nmom_eq; imom++) for (int r=0; r<Nr; r++) s0_ff_eq[iclust][imass][imom][r]=malloc(sizeof(complex)*(3));

        complex *****s0_inv_ff_eq=malloc(sizeof(complex****)*(Nclust+1));
        for (int iclust=0; iclust<Nclust+1; iclust++) s0_inv_ff_eq[iclust]=malloc(sizeof(complex***)*(Nmass));
        for (int iclust=0; iclust<Nclust+1; iclust++) for(int imass=0; imass<Nmass; imass++) s0_inv_ff_eq[iclust][imass]=malloc(sizeof(complex**)*(Nmom_eq));
        for (int iclust=0; iclust<Nclust+1; iclust++) for(int imass=0; imass<Nmass; imass++) for (int imom=0; imom<Nmom_eq; imom++) s0_inv_ff_eq[iclust][imass][imom]=malloc(sizeof(complex*)*(Nr));
        for (int iclust=0; iclust<Nclust+1; iclust++) for(int imass=0; imass<Nmass; imass++) for (int imom=0; imom<Nmom_eq; imom++) for (int r=0; r<Nr; r++) s0_inv_ff_eq[iclust][imass][imom][r]=malloc(sizeof(complex)*(3));


      double *****ZOr_eq=malloc(sizeof(double****)*(5));
      for (int igamma=0; igamma<5 ; igamma++) ZOr_eq[igamma]=malloc(sizeof(double***)*(Nclust+1));
      for (int igamma=0; igamma<5 ; igamma++) for (int iclust=0; iclust<Nclust+1; iclust++) ZOr_eq[igamma][iclust]=malloc(sizeof(double**)*(Nmass_comb));
      for (int igamma=0; igamma<5 ; igamma++) for (int iclust=0; iclust<Nclust+1; iclust++) for(int imass=0; imass<Nmass_comb; imass++) ZOr_eq[igamma][iclust][imass]=malloc(sizeof(double*)*(Nmom_eq));
      for (int igamma=0; igamma<5 ; igamma++) for (int iclust=0; iclust<Nclust+1; iclust++) for(int imass=0; imass<Nmass_comb; imass++) for (int imom=0; imom<Nmom_eq; imom++) ZOr_eq[igamma][iclust][imass][imom]=malloc(sizeof(complex)*(NOa2));

     double ****Zqr_eq=malloc(sizeof(double***)*(Nclust+1));
     for (int iclust=0; iclust<Nclust+1; iclust++) Zqr_eq[iclust]=malloc(sizeof(double**)*(Nmass));
     for (int iclust=0; iclust<Nclust+1; iclust++) for(int imass=0; imass<Nmass; imass++) Zqr_eq[iclust][imass]=malloc(sizeof(double*)*(Nmom_eq));
     for (int iclust=0; iclust<Nclust+1; iclust++) for(int imass=0; imass<Nmass; imass++) for (int imom=0; imom<Nmom_eq; imom++)  Zqr_eq[iclust][imass][imom]=malloc(sizeof(complex)*(NOa2));

     double *aux_mom=malloc(sizeof(double)*Nmom);
     double *aux_eq_mom=malloc(sizeof(double)*Nmom_eq);

    

     for (int sigma=0; sigma<3; sigma++) for (int r=0; r<Nr; r++) for(int imass=0; imass<Nmass; imass++)  for (int iclust=0; iclust<Nclust+1; iclust++) 
     {
		for (int imom=0; imom<Nmom; imom++) aux_mom[imom]=s0_ff[iclust][imass][imom][r][sigma];
		Average_over_Equivalent_Momenta(SinP2,aux_mom,aux_eq_mom);
		for (int imom=0; imom<Nmom_eq; imom++) s0_ff_eq[iclust][imass][imom][r][sigma]=aux_eq_mom[imom];
		
		for (int imom=0; imom<Nmom; imom++) aux_mom[imom]=s0_inv_ff[iclust][imass][imom][r][sigma];
                Average_over_Equivalent_Momenta(SinP2,aux_mom,aux_eq_mom);
                for (int imom=0; imom<Nmom_eq; imom++) s0_inv_ff_eq[iclust][imass][imom][r][sigma]=aux_eq_mom[imom];

     }


  for (int iOa2=0; iOa2<NOa2; iOa2++)
   for (int igamma=0; igamma<5; igamma++)
   {
     int imass=0;
     for (int im2=0; im2<Nmass; im2++) for (int im1=im2; im1<Nmass; im1++)
     {
       for (int iclust=0; iclust<Nclust+1; iclust++)
       {
         for (int imom=0; imom<Nmom; imom++) aux_mom[imom]=ZOr[igamma][iclust][imass][imom][iOa2] ;
         Average_over_Equivalent_Momenta(SinP2,aux_mom,aux_eq_mom);
         for (int imom=0; imom<Nmom_eq; imom++)  ZOr_eq[igamma][iclust][imass][imom][iOa2]=aux_eq_mom[imom];
       }
       imass++;
     }
   }

   for (int iOa2=0; iOa2<NOa2; iOa2++)
   for (int imass=0; imass<Nmass; imass++)
    for (int iclust=0; iclust<Nclust+1; iclust++)
   {
    for (int imom=0; imom<Nmom; imom++) aux_mom[imom]=Zqr[iclust][imass][imom][iOa2];
    Average_over_Equivalent_Momenta(SinP2,aux_mom,aux_eq_mom);
    for (int imom=0; imom<Nmom_eq; imom++)  Zqr_eq[iclust][imass][imom][iOa2]=aux_eq_mom[imom];

   }


     printf("============================\n\n\n");




	/////////////////
	// Write output
	////////////////



     complex **aux_allmomenta=malloc(sizeof(complex*)*(Nmom));
     for (int imom=0; imom<Nmom; imom++) aux_allmomenta[imom]=malloc(sizeof(complex)*(Nclust+1));

     double *aux_clust_real=malloc(sizeof(double)*(Nclust+1));


      
      printf(" WRITING OUTPUT IN %s\n",dir_qav);
      for (int sigma=0; sigma<3; sigma++)
        for (int r=0; r<Nr; r++) for (int imass=0; imass<Nmass; imass++)
        {
        // ff S(p)
	sprintf(file,"ff_r%d_m%d_sigma%d",r,imass,sigma+1);
        for (int imom=0; imom<Nmom_eq; imom++){
          for (int iclust=0; iclust<Nclust+1; iclust++)
                aux_clust_real[iclust]=creal(s0_ff_eq[iclust][imass][imom][r][sigma]);
		write_output_onemomenta_real(dir_qav,file,aux_clust_real,SinP2_eq,imom);
	   }
        sprintf(file,"%s/FF_ALLP/ff_r%d_m%d_sigma%d",dir_qav,r,imass,sigma+1);
        FILE *fout=fopen(file,"w");
        if(fout == (FILE *)NULL) {
              fprintf(stderr,"Unable to open file %s for writing\n",file);
              return EXIT_FAILURE;
         }
         for (int imom=0; imom<Nmom_eq; imom++) {
                for (int iclust=0; iclust<Nclust+1; iclust++) aux_clust_real[iclust]=creal(s0_ff_eq[iclust][imass][imom][r][sigma]);
                fprintf(fout,"%8.10f %8.10f %8.10f  \n",SinP2_eq[imom],creal(s0_ff_eq[Nclust][imass][imom][r][sigma]),error_jack_real(aux_clust_real,Nclust));
          }
      fclose(fout);
     }



      for (int sigma=0; sigma<3; sigma++)
        for (int r=0; r<Nr; r++) for (int imass=0; imass<Nmass; imass++)
        {
        // ff_inv S(p)
        sprintf(file,"ff_inv_r%d_m%d_sigma%d",r,imass,sigma+1);
        for (int imom=0; imom<Nmom_eq; imom++){
          for (int iclust=0; iclust<Nclust+1; iclust++)
                aux_clust_real[iclust]=creal(s0_inv_ff_eq[iclust][imass][imom][r][sigma]);
                write_output_onemomenta_real(dir_qav,file,aux_clust_real,SinP2_eq,imom);
           }
        sprintf(file,"%s/FF_ALLP/ff_inv_r%d_m%d_sigma%d",dir_qav,r,imass,sigma+1);
        FILE *fout=fopen(file,"w");
        if(fout == (FILE *)NULL) {
              fprintf(stderr,"Unable to open file %s for writing\n",file);
              return EXIT_FAILURE;
         }
         for (int imom=0; imom<Nmom_eq; imom++) {
                for (int iclust=0; iclust<Nclust+1; iclust++) aux_clust_real[iclust]=creal(s0_inv_ff_eq[iclust][imass][imom][r][sigma]);
                fprintf(fout,"%8.10f %8.10f %8.10f  \n",SinP2_eq[imom],creal(s0_inv_ff_eq[Nclust][imass][imom][r][sigma]),error_jack_real(aux_clust_real,Nclust));
          }
      fclose(fout);
     }



        // Zq RPROM QAVERAGE
    for (int iOa2=0; iOa2<NOa2; iOa2++)
    for (int imass=0; imass<Nmass; imass++)
     {
        whichOa2(iOa2, Oa2print);
        sprintf(file,"Zq_%s_m%d",Oa2print,imass);
         for (int imom=0; imom<Nmom_eq; imom++){
          for (int iclust=0; iclust<Nclust+1; iclust++)
                aux_clust_real[iclust]=Zqr_eq[iclust][imass][imom][iOa2];
          write_output_onemomenta_real(dir_qav,file,aux_clust_real,SinP2_eq,imom);
         }

         sprintf(file,"%s/Z_ALLP/Zq_%s_m%d.dat",dir_qav,Oa2print,imass);
         FILE *fout=fopen(file,"w");
         if(fout == (FILE *)NULL) {
              fprintf(stderr,"Unable to open file %s for writing\n",file);
              return EXIT_FAILURE;
         }
         for (int imom=0; imom<Nmom_eq; imom++) {
                for (int iclust=0; iclust<Nclust+1; iclust++) aux_clust_real[iclust]=Zqr_eq[iclust][imass][imom][iOa2];
                fprintf(fout,"%8.10f %8.10f %8.10f  \n",SinP2_eq[imom],Zqr_eq[Nclust][imass][imom][iOa2],error_jack_real(aux_clust_real,Nclust));
          }
      fclose(fout);
     }

   // ZO RPROM GAMMAPROM QAVERAGE
   for (int iOa2=0; iOa2<NOa2; iOa2++)
   for (int igamma=0; igamma<5; igamma++)
   {
        whichOa2(iOa2, Oa2print);
        whichGammaProm(igamma,gammaprint);
        int imass=0;
        for (int im2=0; im2<Nmass; im2++) for (int im1=im2; im1<Nmass; im1++)
         {

         sprintf(file,"Z%s_%s_m%d%d",gammaprint,Oa2print,im1,im2);
         for (int imom=0; imom<Nmom_eq; imom++) {
                for (int iclust=0; iclust<Nclust+1; iclust++) aux_clust_real[iclust]=ZOr_eq[igamma][iclust][imass][imom][iOa2];
                 write_output_onemomenta_real(dir_qav,file,aux_clust_real,SinP2_eq,imom);
        }
          sprintf(file,"%s/Z_ALLP/Z%s_%s_m%d%d.dat",dir_qav,gammaprint,Oa2print,im1,im2);
          FILE *fout=fopen(file,"w");
          if(fout == (FILE *)NULL) {
              fprintf(stderr,"Unable to open file %s for writing\n",file);
              return EXIT_FAILURE;
          }
         for (int imom=0; imom<Nmom_eq; imom++) {
                for (int iclust=0; iclust<Nclust+1; iclust++) aux_clust_real[iclust]=ZOr_eq[igamma][iclust][imass][imom][iOa2];
                fprintf(fout,"%8.10f %8.10f %8.10f  \n",SinP2_eq[imom],ZOr_eq[igamma][Nclust][imass][imom][iOa2],error_jack_real(aux_clust_real,Nclust));
        }
      fclose(fout);

        imass++;
          }
   }


     printf("====================================================\n\n\n");

////// FINISH
	printf("Free memory\n");

	free(s0_ff); free(s0_inv_ff);
	free(s0_ff_eq); free(s0_inv_ff_eq);
	free(Zq); free(Lambda); free(ZOr);
	free(Zqr); free(Lambdar); 
	free(Zqr_eq);
	printf("Calc_FF_AmpGF_and_Z.c finished OK ! \n");
	
  return 0;
}
